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ABSTRACT 

O 
O 

Despite observational evidence for cold neutral astrophysical accretion disks, 
the viscous process which may drive the accretion in such systems is not yet un- 
<^ ■ derstood. While molecular viscosity is too small to explain the observed accretion 

efficiencies by more than ten orders of magnitude, the absence of any linear in- 
stability in Keplerian accretion flows is often used to rule out the possibility of 



CN ■ turbulent viscosity. Recently, the fact that some fine tuned disturbances of any 

^ ■ inviscid shear flow can reach arbitrarily large transient growth has been proposed 

as an alternative route to turbulence in these systems. We present an analytic 
psj ■ study of this process for 3D plane wave disturbances of a general rotating shear 

flow in Lagrangian coordinates, and demonstrate that large transient growth is 
the generic feature of non-axisymmetric disturbances with near radial leading 
r£3 ■ wave vectors. The maximum energy growth is slower than quadratic, but faster 

than linear in time. The fastest growth occurs for two dimensional perturbations, 
2 ■ and is only limited by viscosity, and ultimately by the disk vertical thickness. 

After including viscosity and vertical structure, we find that, as a func- 
tion of the Reynolds number, 1Z, the maximum energy growth is approximately 



0.4(7?./ log7?.) 2 / 3 , and put forth a heuristic argument for why 1Z > 10 4 is required 
to sustain turbulence in Keplerian disks. Therefore, assuming that there exists 
a non-linear feedback process to replenish the seeds for transient growth, astro- 
physical accretion disks must be well within the turbulent regime. However, large 
3D numerical simulations running for many orbital times, and/or with fine tuned 
initial conditions, are required to confirm Keplerian hydrodynamic turbulence on 
the computer. 
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1. Introduction 

Accretion disks are very common in astrophysics. Disks are found in Active Galactic 
Nuclei (AGN), around newly forming stars (proto-planetary disks), and surrounding compact 
stellar remnants (white dwarfs, neutron stars, black holes) in binary systems (e.g., Pringle 
1981). Despite the overwhelming evidence for the existence of accretion disks (e.g., Lin & 
Papaloizou 1996), and decades of literature dedicated to their understanding (e.g., see Balbus 
& Hawley 1998, and references therein), the detailed working mechanism of accretion disks 
remains enigmatic at best. 

One of the early puzzles in understanding the accretion phenomenon was the clear 
inadequacy of molecular viscosity in driving accretion in a Keplerian disk. This led to the 
speculation of turbulence as a proxy for viscous transfer of angular momentum (Weizsacker 
1948; Shakura & Sunyaev 1976). The idea was particularly attractive because of the low 
viscosity of astrophysical fluids (high Reynolds number ~ 10 10 — 10 14 ), which, according to 
conventional wisdom, should make the shear flow in accretion disks unstable to turbulence. 

However, in the context of Keplerian disks relevant to most astrophysical applications, 
no instability could be identified. The linear instabilities that seemingly induce the onset of 
turbulence in normal shear flows are stabilized by the Coriolis force associated with rotation 
in a Keplerian disk. Nonlinear hydrodynamic simulations appeared to confirm the absence 
of turbulence in such disks (Balbus, Hawley, & Stone 1996; Hawley, Balbus, & Winters 1999, 
both hereafter referred to as BHSW). 

The breakthrough came in 1991, through re-discovery of the Magnetorotational Insta- 
bility (MRI; Chandrasekhar 1960) by Balbus and Hawley (Balbus & Hawley 1991), who 
showed through Magnetohydrodynamic (MHD) simulations that initial seed magnetic fields 
in an MHD disk grow exponentially within a few rotation times, causing the onset of MHD 
turbulence. The MRI instability is now widely accepted as the driver of MHD turbulence in 
ionized disks (see e.g., Balbus 2003, for a recent review). 

Despite the great success of the MRI in many astrophysical systems, it is known that this 
instability will not operate in disks with very small ionization fractions, where the magnetic 
flux is poorly coupled to the gas. Examples of systems with low ionization fractions are 
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proto-planetary disks, outer regions of AGN disks, and white dwarf disks in the low state 
(Gammie 1996; Gammie & Menou 1998; Formang, Terquem, & Balbus 2002). The route to 
turbulence and accretion in such neutral disks remains an outstanding puzzle in theoretical 
astrophysics. 

On the other hand, laboratory experiments of Taylor- Couette systems (fluid flow be- 
tween concentric rotating cylinders) seem to indicate that, although Coriolis force delays the 
onset of turbulence, the flow is ultimately unstable to turbulence for Reynolds numbers larger 
than a few thousand (e.g., see Richard 2001), even for subcritical systems (systems with no 
linear instability). Longaretti (2002) reviews the experimental evidence for the existence of 
turbulence in subcritical laboratory systems, and based on phenomenological analogy, con- 
cludes that a similar process must happen in astrophysical accretion flows. Longaretti (2002) 
also claims that the absence of turbulence in previous numerical simulations (BHSW) is due 
to their small effective Reynolds number, which is limited by the numerical viscosity caused 
by the finite resolution of the simulation. Indeed, Bech & Andersson (1997) see turbulence 
persisting in numerical simulations of subcritical rotating flows for large enough Reynolds 
numbers. 

How does a shearing flow that is linearly stable to perturbations switch to a turbu- 
lent state? A possible explanation, known as bypass transition, has been discussed in the 
fluid mechanics community for some time (see Grossmann 2000; Reshotko 2001; Schmid & 
Henningson 2000, and references therein, for an overview), though its diffusion into the as- 
trophysical community has been slow (Ioannou & Kakouris 2001; Chagelishvili et al. 2003; 
Tevzadze et al. 2003; Yecko 2004; Umurhan & Regev 2004). The bypass concept is based on 
the fact that definite frequency linear modes are not orthogonal in a shear flow. Therefore, 
even if all of the linear modes are decaying, a suitably tuned linear combination of them can 
still show an arbitrarily large transient energy growth in the absence of viscosity. In lieu of 
linear instabilities such as MRI, the transient energy growth, supplemented by a non-linear 
feedback process to repopulate the growing disturbances, could plausibly sustain turbulence 
for large enough Reynolds numbers. 

This paper, along with a companion paper (Mukhopadhyay, Afshordi, & Narayan 2004, 
hereafter MAN04), investigates the transient growth of perturbations in rotating shear flows, 
with an emphasis on applications to astrophysical accretion disks. Both papers study the 
dependence of the transient energy growth on various parameters of the system, in particular 
the wave vector of the perturbations and the Reynolds number. While MAN04 focuses 
on an eigenmode analysis in Eulerian coordinates for a shearing flow restricted between 
rigid walls, the present paper studies the bypass process in Lagrangian coordinates for an 
infinite shearing flow. The advantage of the Lagrangian approach is that there is no explicit 
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coordinate dependence in the equations, and therefore the solution can be decomposed into 
plane waves. The drawback, however, is the explicit time dependence of the equations which 
prohibits definite frequency solutions, and thus requires explicit integration in time for each 
mode. The two analyses presented here and in MAN04 involve different numerical/analytic 
techniques and are both, in our opinion, valuable. The consistency of the results validates 
the general picture of our understanding of transient growth in the linear regime. 

In §2 we summarize the current understanding of the transient growth phenomenon 
in shearing flows. §3 introduces the basic linear equations in Lagrangian coordinates in 
their most general form. §4 constitutes the main body of the paper, where we approach 
the problem of transient growth in inviscid and incompressible flows for general plane wave 
solutions. This is followed by §5 and §6, where we study the effects of viscosity and compress- 
ibility/vertical structure, respectively, and make the connection to astrophysical accretion 
flows. Finally, in §7 we discuss conditions for the emergence of turbulence, and its realization 
in numerical simulations. §8 summarizes our results and concludes the paper. 



2. Understanding Transient Growth through Swinging Plane Waves 

As mentioned in §1, the systematic approach to transient growth in subcritical systems is 
through a linear combination of non-normal decaying modes. In particular, in a local region 
of an accretion flow, different definite frequency modes with equal vertical and azimuthal 
wave numbers, but different radial profiles, are generically not orthogonal to one another, 
i.e. 

J d 3 x 5v a • 5v b ^ 0, (1) 

where <5v a and 5v b are velocity profiles of different modes. Therefore, even though all the 
modes may decay with time, a solution may still show a temporary energy growth for suit- 
able initial conditions because of the cross terms in the energy expression. For a detailed 
description of the eigenmode approach, we refer the reader to MAN04 and Yecko (2004), 
which explain the numerical methods used to find the maximum growth through optimizing 
a linear combination of eigenmodes. 

An alternative approach, which allows simple analytic treatment of linear perturbations 
for the case of astrophysical accretion flows, is the so-called shearing box approximation in 
which we study the linear evolution of plane wave perturbations in Lagrangian coordinates 
within a small region of an accretion flow (Chagelishvili et al. 2003; Tevzadze et al. 2003; 
Umurhan & Regev 2004). This is the approach we follow in this paper. 

Fig. 1 shows our choice of local Cartesian coordinates: x is along the radial direction, 
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y is along the azimuthal or streamwise direction, and z is along the vertical direction. The 
unperturbed flow has a velocity in the y-direction and a velocity gradient (shear) along the 
x direction. The Coriolis force associated with rotation is described by an angular frequency 
vector f2 pointed in the z direction. We define 

q = -dlnSl/dlnR, (2) 

which is a dimensionless parameter that quantifies the shear in the local comoving box 1 . 
For example, q = 3/2 corresponds to a Keplerian accretion disk and q = 2 to a disk with 
constant specific angular momentum. 

Starting with a plane wave in Lagrangian coordinates with radial wave number k% (Eq. 
16), we show in §3 that the Eulerian radial wave number k x of the plane wave evolves as a 
function of time t according to: 

K = k^ + { q nt)k y , (3) 

while the azimuthal and vertical wave numbers k y and k z remain unchanged. Thus, the 
plane wave is effectively frozen into the flow and is swung around by the shear. For sim- 
plicity, as in Chagelishvili et al. (2003) and Umurhan & Regev (2004), let us consider 'two- 
dimensional' (k z = 0) incompressible and inviscid perturbations. In this regime, in general, 
the two-dimensional vorticity £ = k x 5v y — k y 5v x remains constant, while the two-dimensional 
divergence vanishes due to incompressibility, i.e. k x 5v x + k y 5v y = 0. Therefore, the energy 
content of the plane wave scales as: 

^^-^-rc + J^' + s- (4) 

which reaches a maximum when k x = 0, or k^/k y = —qVLt. If the initial wave vector k x = k% 
is negative, then the maximum is reached at positive time, and the energy growth factor is 
given by: 

G max = |^ ~ (k L Jk y f = (qnt) 2 for - k L Jk y » 1. (5) 

As promised, the growth can become arbitrarily large for long enough time. As we show 
in §6, the maximum growth is limited only by viscosity and scales as 7£ 2//3 , where 1Z is the 
Reynolds number (Chagelishvili et al. 2003; Yecko 2004, MAN04). Note that, after reaching 
the maximum, the linear solution drops and decays to zero asymptotically. Therefore, to 
have sustained turbulence, the perturbations must become non-linear before reaching the 
peak and must provide sufficient feedback to keep the perturbations going. Umurhan & 
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Regev (2004) describe a 2D numerical simulation of the local accretion flow in which they 
find sustained turbulent behavior in the absence of viscosity. 

Moving away from the k z — plane, it is seen that for q < 2 (the regime of interest 
for astrophysical disks), the maximum growth drops (Yecko 2004, MAN04). However, 
Tevzadze et al. (2003) argue that for k z ^ 0, although the maximum growth is smaller, 
vertical stratification may cause the solution to have a non-vanishing asymptotic value, 
which is a fraction of the maximum growth. It is not clear if this will significantly help the 
onset of turbulence. 

This concludes our summary of (the rather thin) astrophysical literature on transient 
growth and the bypass mechanism. In the following sections, we present a formal treatment 
of the Lagrangian hydrodynamic equations in the shearing box approximation, and study 
the transient growth of general plane wave solutions. 



3. Hydrodynamic Equations in the Local Lagrangian Coordinates (Shearing 

Box Approximation) 

We limit the calculations to scales much smaller than the thickness of the disk, which, 
for a geometrically thin disk, is in turn much smaller than the distance to the central object. 
We thus ignore the boundaries. Within a comoving local box, and in terms of Lagrangian 
coordinates, the Navier-Stokes and continuity equations can be written as: 

v = -c 2 s VA + zA7 2 v + 2vxH, (6) 

A = -V.v, (7) 

r = v(r L ), (8) 
c)r L 

V S "L.V 1 , (9) 

where v is the fluid velocity vector, r and r L are Eulerian and Lagrangian position vectors 
respectively, v is the kinematic coefficient of (molecular) viscosity, and A = lnp is the 
logarithm of fluid density. Since the local box rotates with the flow, we have a Coriolis term 
on the right hand side of the Navier-Stokes equation (6), which is proportional to the local 
angular velocity Q; in the shearing box approximation, the Coriolis parameter Q is taken to 
be independent of position. We have dropped the centrifugal and gravitational accelerations, 
as they cancel in the equilibrium flow and do not contribute to the perturbation equations. 

Fig. 1 shows the unperturbed velocity field within the local comoving box. Let us define 
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f2q as minus the gradient of the unperturbed velocity field, 

o q n o \ 

Oq = -Vv: HW)xR=| ); g = (10) 







where R = (i?, 0,0), and i? is the distance to the center of the disk. Notice that, as the 
velocity field is normal to its gradient, and flow lines are straight in the local approximation, 
Eq. (8) can be integrated to give 

dr L 

r L = r + ntr.q — = 1 + tttq, (11) 
or 

and thus 

V = (l + fitq).V L . (12) 

In deriving the equations for Eulerian perturbations in Lagrangian coordinates, we note 
that, as the unperturbed velocity field v has a spatial gradient (Fig. 1), it has a non- 
vanishing time derivative in the perturbed Lagrangian coordinates, i.e. 

5\ = v — v = v — v.Vv = v + f2v.q. (13) 

Combining this with Eqs. (6) and (7) yields the linear perturbation equations: 

<5v = -c 2 s X75\ + z/V 2 5v + 25v x Q + IWv.q, (14) 
SX = -V.5v, (15) 

where the Eulerian gradients are related to the Lagrangian gradients according to Eq. (12). 
We note that since q is constant, the Navier-Stokes and continuity equations have no explicit 
dependence on the Lagrangian coordinates. Therefore, we can decompose a general linear 
perturbation into plane wave solutions of the form: 

5v, 5\ocexp(ik L .r L ). (16) 

Note that, although the Lagrangian fluid equations (6 and 7) look linear in v and A, the 
above plane wave description breaks down for non-linear perturbations, since Eq. (12) will 
be modified at higher orders. 

For a Lagrangian plane wave solution, the Eulerian wave number k can be obtained 
from Eq. (12): 

k = (k x , ky, k z ) = (1 + fitq).k L = {k L x + qntk L y , ky , k L z ). (17) 



- 8- 



This can be combined with Eqs. (14) and (15) to yield the mode equations for Lagrangian 
plane wave solutions 

5v x = -ic 2 s (k^ + Qtqk^)5X + 2Q5v y - uk 2 5v x , (18) 
5v y = -ic 2 s ky5X + (q- 2)Q5v x - uk 2 5v y , (19) 
5v z = -ic 2 s k^5X-uk 2 5v z , (20) 
5X = -i\{k L x +ntqk^)5v x + k^5v y + k^5v z }, (21) 

where k is the total Eulerian wave-vector 

k 2 = {k L x + ntqk L y f + {k L y f + {k L z f. (22) 

For given values of q, f2, z/, and c s , Eqs. (18-21) provide a homogeneous set of linear first 
order differential equations with time dependent coefficients. In the rest of this paper, we 
attempt to study the growth of energy £, defined as 

£ = |(KI + KI + M), (23) 

for different modes in various limits of the parameter space of interest to astrophysical/physical 
problems. As we discussed in §2, the presence of a large transient energy growth, even in 
a small region of phase space, may act as a possible trigger for the onset of self-sustained 
turbulence. 

Eqs. (18-21) can be simplified by introducing the 2D divergence and vorticity, A and 
£, defined as 

A = k x 5v x + k y 5v y , (24) 
£ = k x 5v y - k y 5v x , (25) 

which yield 

i = (q-2)SlA-vk% (26) 

A = (w^-" k2 ) & + 2n ( 1 -WT^) i - ic2 - ik ' + k ' )sx ' (27) 

\ x jj / v x y / 



y 

5v z = —ic 2 k z 5X — uk 2 5v z , (28) 

5X = -i[A + k z 8v z ]. (29) 

In terms of the new variables, S can be written as 

£ -l(kl + W + 4 ' (30) 
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4. Idealized Inviscid and Incompressible Flow 



The Reynlods number corresponding to molecular viscosity is typically very large in 
astrophysical accretion disks. In this limit, we can visualize a regime in which 



(c s t)- 2 < k 2 < iyt)~ x , 



(31) 



where t is the characteristic time of the transient growth. Within this regime, we may neglect 
viscous effects [y ~ 0) and we may also consider the fluid to be effectively incompressible, 
i.e. 5\ — > and c 2 s — > oo, while the pressure perturbation c 2 5\ remains finite. 

We begin by rearranging Eqs. (26-29) to get: 



i = (q-2)nA-uk% 



A = 



k 2 

A 4- I — 

k\kl + k 2 ) \ k 2 



qk 2 



(k 2 x + hi) 



k 2 z (ic 2 5\) = A + uk 2 A-i(5\ + uk 2 5\), 
k z 5v z = —A + i5X, 



i-uk 2 A + l [ k ^pL 



(32) 

{5X + uk 2 5X), 

(33) 
(34) 
(35) 



with the energy given by 



c _ |A 2 | + |e 2 | |A| 2 + |5'A| 2 



m + k 2 ) 



Ak 2 



(36) 



In the incompressible and inviscid limit, we may ignore the terms linear in 5 A and v on 
the right hand sides of Eqs. (32-36). Eqs. (32) and (33) then take the form 



i = (q-2)nA, 

: 2ttk 2 

A = 



k 2 (k 2 x + k 2 ) 



{qk x k y A + [{l-q)k 2 y + k 2 x ]i}. 



(37) 
(38) 



As the coefficients in the above equations are real, we can without loss of generality focus 
on real solutions, and thus the energy becomes 



2 cl 



£ = 



k 2 A 2 + k 2 C 



(39) 



where k x = k x + (qflt)k y , is the only time- dependent component of the wave vector. 
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4.1. Near two-dimensional Growth (k z <C k y ,k x ) 



Eqs. (37-39) are significantly simplified in the two-dimensional regime, k z = 5v z = 
(see Chagilishvili et al. 2003; Umurhan & Regev 2004). In this limit, we find 



A = 

£-2D = 



5X = 0, £ = const., 
£ 2 £ 2 



(40) 
(41) 



We have already seen in §2 that this leads to a maximum energy growth factor of (qVlt) 2 for 
qVLt = -k L x jk y > 1 (Eq. 5). 

Starting with A = and £ = £ = const, as the zeroth order solution, to first order in 
k 2 , Eq. (38) becomes 



2ttk 2 z 

U2 _i_ 1.2 
tv x -|- tv y 



1 - 



qkl 



V2 I 1.2 



£o, 



which can be integrated to give 



A = A + 



U2 

K y 



{2q~ 1 - 1) tan" 1 (k x /ky) 



kry-k 



kl + k 2 y 



+ 0(k 4 z ). 



(42) 



(43) 



Here A = 0{k 2 z ) and so it vanishes in the limit k z — > 0. Plugging this back into Eq. (37), 
yields £ to first order in k 2 



£ = £o + (<?-2)fiA t--||(2-<?) 



A- 



i(k x . 



( 2 -g)^tan- 1 (-^)-ln 1 + -f 



A: 



A'v, 



A-: 



?,2 



(44) 



Now we can go back to Eq. (39) to find the evolution of energy with time. Focusing on the 
solution with k z = and maximum growth with k x = —(qQt)k y , after some manipulations 
and marginalizing over £ and A , we end up with an expression for the modified maximum 
growth 



G max (t) = ( q ntf 



1-4 



qk y 



(2 - q) Infant) + O {k 2 Jkl) 



(45) 



We see that, as expected (see §2), the growth function decreases as we move away from the 
k y axis so long as q < 2. 



4.2. Near Axi-Symmetric Growth (k y <^ k Z) k x ) 



Another special case which allows a simple analytic solution is the case of axi-symmetric 
perturbations, where k y = 0. Since there is no azimuthal dependence, the wave pattern is not 
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swung round by the unperturbed flow, and therefore, there is no explicit time dependence 
in the equations (see Eq. 17). As a result, Eqs. (32-36) simplify to 

i = -(2-q)QA, (46) 

= W^' (47) 

which allow a harmonic solution 

A = Aexp(-iujt) + Bexp(iut), (48) 

f = i — ~ 2 ^ Q [Aexp(-iujt) - B exp(iut)], (49) 



where 



_ 2(2 -q)tfkl 

" kl + kl ■ (50) 



We first note that for q > 2 the frequency uo = ±i\u>\ is imaginary. This means that 
the system has exponentially growing perturbations, a reflection of the Rayleigh stability 
criterion (according to which a flow with specific angular momentum decreasing outward is 
unstable). However, this regime is not of interest since real disks cannot survive here. 

Let us, therefore, consider q < 2, where stable circular orbits are allowed and an accre- 
tion flow is possible. Without loss of generality, we can change the origin of time to eliminate 
the phase difference between A and B. The leftover phase will be irrelevant for calculating 
the energy, and thus we can assume that A and B are real. Plugging this into Eq. (36), 
after straightforward manipulations, we arrive at 

S AS = 1 -{K 2 + k; 2 ) [(A 2 + B 2 ) (4 - q) + 2qABcos(2out)} , (51) 

where the subscript AS identifies the axi-symmetric solutions. 

Since S is a periodic function of time, the maximum growth factor is equal to the ratio 
of its maximum to its minimum: 

_ (A 2 + B 2 ) (A-q) + 2AB\q\ 
GAS ~ (A^TW) (4 - q) - 2AB \q\' [ } 



which is maximized for A = B: 



sgn(g) 



GAS,max = ~. i 7 = ~ • (53) 

A-q-\q\ \2-qJ 

For most astrophysical disks, q > 0, i.e. the angular velocity decreases with increasing 
radius, and therefore significant growth can only happen in the limit of a constant angular 
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momentum disk (2 — q = d\n(HR 2 )/d\nR — > 0). For a Keplerian disk with q = 3/2, the 
maximum possible energy growth is only a factor of 4. 

Focusing on the marginal case of a disk with constant angular momentum (g = 2), since 
the harmonic frequency u is zero, we find solutions that are linear in time: 

2flk 2 z & 



A = 



£-2 _i_ U2 ' 



£ = const.. 



The growth is then given by 



G AS {t;q = 2) = 1 + 



-(2fit) 5 



(54) 



(55) 



L2 _i_ J.2 

As in the case of two-dimensional perturbations, we see that the growth for axisymmetric 
perturbations of a constant angular momentum disk is quadratic in time and can be arbi- 
trarily large in the absence of viscosity. In contrast, for any q < 2, there is a firm limit to 
the maximum growth, given by equation (53). 

Now, let us consider small but non-vanishing values of k y . As a small value of k y 
corresponds to a slowly changing frequency, we can use the WKB approximation to solve 
Eqs. (37-38). This gives 



£(t) ~|(i)exp 



dt u(t) 



, and A(f) ~ A(t)exp 



dt u(t) 



Ignoring terms of order k 2 in Eq. (38), we find after some simple manipulations that 

d 90P 
A = A|ln(Jb^) + ^. 

Plugging (56) into Eqs. (37) and (57), and keeping only terms linear in k y , yields 

kl 



iu(t)£ ~ (2 - q)£lA , and £ 2 oc 



cuk 



2 ' 



(56) 



(57) 



(58) 



Plugging this result into the expression for energy (Eq. 39), and ignoring second order terms 
in k y , we find 

£wkb oc k~ l = {k 2 x + kl)- 1 / 2 , (59) 
which reaches a maximum when k x = k% + k y {qVti) = 0. Therefore, the maximum growth 
becomes 



m 2 + ki 



~ (qnt) 



1 + 



for qQt > 1, 



(60) 



Thus, even near axi-symmetric plane-wave solutions can achieve arbitrarily large tran- 
sient growth. However, the growth is only linear in time, whereas it is quadratic for 
two-dimensional solutions. Therefore, the growth is significantly slower than for the two- 
dimensional case. 
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4.3. Ideal Growth for a General Plane Wave Solution 

Having analysed the regions near the two axes in wave vector space, we now consider 
general perturbations of an incompressible and inviscid flow. We have not been able to find 
an analytic solution to Eqs. (37-38), but we are able to put a lower bound on the growth. 

First, consider how the combination £iA2 — £2^1 evolves with time, where (£1, Ai) and 
(£2, A2) are two independent solutions of Eqs. (37-38). From these equations we obtain 

which can be easily integrated to give 

(£1 A 2 - £ 2 A0 (7^2) = const. (62) 

Now consider two solutions with initial conditions (£ , — A ) and (£ ,A ), which have 
the same initial energies: 

where 

kl = (k L x f + k" y + k\. (64) 
After time t, the sum of the energies of the two solutions is: 

k\A\ + A*) + kHil + £1) (**A| + m) + (k*Al + m) 



S 1 (t)+S 2 (t) = 



Aki{ki+kD ±mi+kD 



> 2fc 4^ + |) Al) = l±kl{kl - kl)] 1 (Jkgjfc./Jfc), (65) 

where we have used Eq. (62) to substitute for £iA 2 — £2Ai in terms of its initial condition. 

We can now translate this result to a lower bound on the average growth function of 
the two solutions: 

Gi(t) + G 2 (t) _ S 1 (t) + £ 2 (t) ( 2k z kfe A /k \ 

— I 7.2 A 2 , 7.2 2 • V* 3 "/ 



2 £ 1 (0)+£ 2 (0) ~ \ klAl + klH 

Maximizing this lower bound over the initial conditions (£0, A ), and noting that the average 
of Gi and G 2 is a lower bound on the maximum growth, we finally arrive at: 

G max (t) > Max^ ~ ( q nt)-^==, for qQt » 1, ^. (67) 
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The maximum is achieved for k x ~ —{qVLt)k y . 

We first note that the above lower bound is smaller than the maximum growth near the 
k y axis (Eq. 5), but becomes asymptotically equal to the growth close to the k z axis (Eq. 60). 
More generally, this result implies that arbitrarily large transient growth is a generic feature 
of plane waves with relatively large leading radial wave numbers (k% ^> k y , k z \ k^k y < 0). 

One interesting observation is how the slope of the energy growth function depends on 
position in the k y — k z plane. We find: 

d\nG max (t) \ 1 if k y < k z , 



dln(qttt) \ 2-4q- 1 (2q- 1 -l)(k 2 z /k 2 y ) if k z < k y , ^ 

which is a result of Eqs. (60) and (45). We see that, in general, the growth function behaves 
as a power law in qflt, i.e. G oc (qQt) a , where 1 < a < 2. We can justify this conjecture by 
plugging a scaling ansatz into Eqs. (37-38): 

$ = |A£, A = Akl (69) 
This ansatz satisfies the equations in the limit k x ^> k y ,k z , if 

P = 7 + 1, (70) 
7(7+1) = -2q-\2q- 1 -l)(k z /k y ) 2 , (71) 



or 



2 7 = -1 ± yJl-Sq-\2q-^-l){k z /k y f. (72) 

Plugging this result into Eq. (39), we see that £ oc k^ 2 ^ ■ Since the scaling solution breaks 
down when k x ~ k y , k z , the maximum growth relative to the initial value of k^ = —(qQt)k y 



takes the form 

G m ax if) — 

where 



{qnt) 1 7PTF 1 F{W 



for iqQt) > 1 (73) 



a = Ma, Re(- 27 ) ^{ \ + ^ _ « £ > £• ,74) 

where 

fx = q [8(2 - q)Y 1/2 ~ 0.20 for g = 3/2. (75) 



We note that the above result for the logarithmic slope of the maximum growth function 
is consistent with our previous asymptotic results for the two limits k y <C k z and k z <C k y 
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(Eq. 68). The exact value of the dimensionless factor F(k z /k y ) ~ 1 does not come out of the 
scaling argument. However, based on our solution in the large and small k y /k z limits (Eqs. 
60 and 45), we can say that F goes to 1 as its argument goes to zero or infinity. Moreover, 
the lower bound in Eq. (67) requires that F(x) > 1 for x > /i. 

Fig. 2 shows how the logarithmic slope of the growth function, a, depends on k z /k y for a 
Keplerian accretion flow. For k z < fik y , the slope monotonically increases with decreasing k z , 
reaching its maximum when k z = 0, as pointed out by Tevzadze et al. (2003); Yecko (2004); 
and also MAN04. Thus, the fastest growth is achieved for two-dimensional solutions. For 
k z > fik y , as for the k z ^> k y case (§4.2), the solution is oscillatory with a growing amplitude, 
yielding a maximum energy growth which is linear in time. 



5. Dependence of the Maximum Growth on Viscosity 



In this section, we study the effect of viscosity on the growth of incompressible local 
perturbations of an accretion disk. Neglecting the terms of order 5X in Eqs. (32-33), but 
keeping the terms that include viscosity, we end up with 



A + uk 2 A 



(q-2)QA, 
2Qk 2 z 



{qk x k y A + [(1 - q)k 2 y + kl]i) . 



(76) 
(77) 



Now, introducing 



i = f • exp 



v / k 2 dt 



, and A = A • exp 



/ 



v \ k 2 dt 



(78) 



we see immediately that £ and A satisfy the same equations that £ and A satisfy (Eqs. 
37-38) for inviscid perturbations. As the expression for the energy (39) remains unchanged, 
we see that in the presence of viscosity, the growth function is simply modified to 



G max {t) = G max (t; u = 0) exp 



t 

' I k 2 dt' 
'o 



(79) 



For qQT ^> 1, based on the results of the last section, we know that the inviscid 
maximum growth is proportional to (qQt) a , where 1 < a < 2, and the maximum is achieved 
when k x = k x + (qQt)k y ~ 0. Plugging these into the above expression for the modified 
growth yields 

2uk 2 

G ma x(t) oc (gfit) CT exp 



3q 



-(qQty 



(80) 
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which reaches a maximum for: 

= V3 • (81) 

Plugging this back into Eq. (79), we find that 

Gmax = G max (t m ax] V = 0) • exp(-(j/3), (82) 

where G max (t max ; v — 0) was obtained in Eq. (73). 

Consider as an example the case of two-dimensional perturbations, where the solution 
has a closed form and the inviscid growth is equal to (cflt) 2 . In the presence of viscosity we 
find 

G max (k g = 0)~(jj^j exp (-2/3) . (83) 



6. Vertical Structure and Finite Speed of Sound 

The heuristic argument for assuming incompressibility in the analysis so far is that, 
if the wavelength of the velocity perturbations is much shorter than the sound horizon for 
the time of interest, then the density perturbations (i.e. sound waves) reach equilibrium 
early on and thus the density is effectively uniform during the timescale of interest for 
velocity perturbations. For a geometrically thin disk around a gravitating mass, the vertical 
half thickness of the disk H is comparable to the sound horizon corresponding to one disk 
rotation time (e.g., Pringle 1981): 

H ~ c s tt-\ (84) 

Therefore, for processes that take longer than one rotation time, wavelengths shorter than 
the disk thickness can be approximately treated as incompressible. 

We can refine this heuristic picture by focusing on the two-dimensional perturbations 
(k z = 0) discussed in §4.1, since their solutions are available in closed form, and solving 
Eqs. (26-29) to first order in c~ 2 . As we have already discussed the effect of viscosity in 
the last section, we will for simplicity ignore viscosity in the following analysis. Assuming 
kz = 8v z — v — 0, Eqs. (26-29) can be simply combined to give 

A = -p^ir < 85 > 

{ = (^ 1 -'* J )e + 2fl'( t -2)(l-^ ? ){ + «*S + ^)(6-{), (87) 

\ x y / \ x y / 
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where £ is a constant. 



Recognizing that £ = £ an d A = <5A = are the solutions to zeroth order in c s 2 , we 
can easily write down the first order solutions: 



= 

A = 



1 + 



2Q 2 (2-q) 



2ifl£ 



C 2 s (kl + k 2 ) 

4qn 2 k x ky£ 
c 2 (k 2 x + k 2 y f 



1 - 



1 - 



qk 2 



qk 2 



k 2 x + k 2 



Co, 



V.2 I U2 

K x -|- fi, y 



[(2? - - m 



(88) 
(89) 
(90) 



Noting that the maximum growth occurs when k x = 0, we find the correction to the growth 
function: 



Gmax (t) 



where H = c s /Q. 



(qVttf {1 - 2(2 - q){q - l)(c s k y /Q)- 2 + 0(c s k y /Q)- 4 } 
(qnt) 2 exp {-2(2 -q)(q- l)(k y H)- 2 + 0(k y H)- A } , 



(91) 



At this point, we should note that the vertical structure of the disk puts a lower limit 
on the vertical component of the wave vector k z ^ min ~ n/(2H) (assuming free boundary 
conditions at z — ±H) 2 . Therefore the pre-factor of (qVlt) 2 must be replaced by (qVtt) a ', 
where a is given in Eq. (74), and is smaller than 2 for a finite value of k z . Let us also define 
the Reynolds number 1Z by 

(92) 



K = 



v 



Modifying Eq. (91) for k z = n/(2H), and using Eq. (81) we arrive at: 

\nG rnax (k z = ~2jj) — 

2 -[l - 2(2 - q)q-\K/2f (k y H)- 2 ] In [qTZ^H)- 2 } - 2 - - 2(2 - q){q - l)(k y H)- 2 + 0(k y H)-\ 

(93) 



2 Since the z-dependent parameter that appears in the equations is c s , H is the characteristic scale for 
variations in c s , and not the un-perturbed density. In particular, for an isothermal disk (c s = const.), k z has 
no minimum, even for a disk with finite thickness. However, for a realistic geometrically thin and optically 
thick disk, c s is expected to drop significantly at the disk surface, which puts a lower limit on the vertical 
wave number k z . 
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Due to the non-algebraic dependence of G max on k y , one cannot express the value of k y 
that maximizes G max , or G max itself, in a closed form. However, we can find an asymptotic 
expansion, in the limit of In ft ^> 1. In this limit, we find 



G ~ 



2exp(— l)g 3 
(2-g)vr 2 J 
ft 



2/3 



n 

toft 



2/3 



0.36 



2/3 



• exp 



tog ft 

This maximum is achieved for 

7r 2 (2-g) 



o 



■ exp 

In In ft 
In ft 



In In ft 
In ft 



, for a Keplerian disk (q — 3/2). (94) 



2q 2 



toft + 0(1) 



~ 2.51ogft + 0(l), for a Keplerian Disk (g = 3/2). 



(95) 



For example, to reach a maximum growth of ~ 10 3 , we need a Reynolds number ~ 10 6 , 
and the maximum growing perturbations have k y H ~ 4, k z H ~ ir/2. We note that, in this 
regime, ignoring the terms of order (k y H)~ 4 , as we did in Eqs. (91-93), introduces < 1% 
error, and is therefore justified. 



7. Discussion 

Throughout this paper, we have demonstrated that arbitrarily large transient growth 
of plane wave perturbations in a stable cold disk is possible in the restricted phase space of 
modes with relatively large radial leading wave numbers (i.e. k x ^> k y ,k z ;k x k y < 0). This 
growth is only limited by viscosity, and eventually by the vertical thickness of the accretion 
disk, as we saw in §5 and §6. 

So why do 3D simulations of local hydrodynamic accretion flows fail to see the onset 
of turbulence for a Keplerian disk (BHSW)? As we mentioned in §1, one possibility is the 
limited numerical resolution of the 3D hydrodynamic simulations (e.g., Longaretti 2002). 
Balbus (2004) invokes the scale invariance symmetry of inviscid Navier-Stokes equations, to 
argue that the lack of any instability on large scales probed by simulations implies stability 
for smaller scales. Although this argument may hold for positive eigenvalue instabilities such 
as MRI, the amplitude of transient growth is directly affected by the dynamical range of the 
simulation (i.e. the effective Reynolds number). Therefore, the resolution of a simulation 
may be a critical factor which decides if an energy growth large enough to sustain turbulence 
can, or cannot be achieved. 
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A yet more important factor in simulating the bypass phenomenon may be the assumed 
initial conditions of the simulations. Although the steady turbulent phase is expected to be 
independent of initial conditions, due to the restricted nature of modes with large transient 
growth, the time needed to reach the steady turbulent phase will significantly depend on 
the choice of initial conditions. For example, based on the results of §4.3, we can show that, 
starting with an isotropic distribution of energy in (3D) k-space, the total linear energy 
decays in the linear regime. 

In order to see this, let us consider the general case of a D-dimensional isotropic initial 
energy distribution within the phase space of modes, i.e. dS = /(A;)d D k. Here, D = 3 
corresponds to a 3D isotropic initial condition, D = 2 is a 2D isotropic distribution with 
k z = 0, and D = 1 refers to initial conditions with k y , k z <C k x . Since, the Lagrangian modes 
are orthogonal, the maximum energy growth is the sum of the maximum energy growth of 
the individual modes, i.e. 

_ fG max (k)f(k)d°k 
Umax ~ Jf(k)d°k ■ W 

At time t ^> the maximum growth is peaked at — k y = and k z = with 

G max ~ (qQt) 2 . Based on the analyses of §4.3, we can also roughly estimate the width of 
the peak, i.e. Ak y ~ (qVLt)~ l k y ~ (qQt)~ 2 k^ and Ak z ~ k y ~ (qflt)~ 1 k^. In particular, for 
isotropic 3D initial conditions, we find 

< w jimwrn (qat) • (97) 

implying that, despite the presence of a large transient growth oc t 2 , as a result of the 
shrinking phase-space volume (which decays as Ak y Ak z oc t~ 3 ) the total energy in the 
perturbations decays as which is completely consistent with the BHSW results. 

In 2D, as it was recently pointed out by Johnson & Gammie (2005), the shrinking of the 
phase space volume exactly cancels the transient growth, which implies a constant energy 
and thus, lack of any significant growth in the linear regime: 

, S(qntf;(ki)( q m)--'k L ,dki _ 

* ,2D ~ — Timwi • ( ' 

which is also consistent with simulations of Umurhan & Regev (2004). 

Therefore, it is only with near one-dimensional initial conditions, i.e. k^ ^> k y , k z , where 
we can obtain significant transient growth of G max ^D = (qflt) 2 until k^ = —(qflt)k y . 

One may speculate that in a 3D simulation the total energy decays until/unless it is re- 
distributed into near radial leading modes through non-linear couplings. BHSW follow the 
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evolution for only a few orbital times, which is probably not enough to reach the expected 
turbulent phase with near-radial structure. However, a direct way to test the viability of 
this scenario is to study a 3D simulation with near-radial initial conditions, and see if a 
self-sustained turbulent phase can be realized. 

A more subtle question to address is the minimum energy growth (or Reynolds number) 
required for the onset of sustained turbulence. Unlike linear evolution, the answer to this 
question requires solving non-linear equations, which can be usually done only through nu- 
merical methods. For example, Umurhan & Regev (2004) see sustained turbulence through- 
out the duration of their inviscid 2D simulation, while their run with a 1Z ~ 10 5 decays in 
a few hundred rotation times. However, it can be formally shown that due to the steady 
decay of comoving vorticity in 2D, one cannot sustain turbulence in a periodic box with- 
out external forcing. Therefore 2D simulations will not be able to pinpoint 1Z C , the critical 
Reynolds number necessary for the onset of sustained turbulence. However, this does not 
necessarily pose a problem for real disks, as the modes that show the largest growth are 
inherently three-dimensional (see the end of §6). 

Assuming that there exists a non-linear feedback process to repopulate the growing 
disturbances, let us present a heuristic way for estimating 1Z C . Most theoretical studies of 
2D turbulence are based on describing the turbulent flow as a gas of 2D interacting vortices 
(see e.g., Miller 1990, and refrences therein). While in isotropic turbulence the vortices are 
typically circular, as can be seen in simulations such as that of Umurhan & Regev (2004), in 
a shear flow, the vortices are stretched along the stream (our y direction). The aspect ratio 
of vortices in their simulation of a 2D Keplerian flow is within the range 6 — 11. Let us now 
conjecture that this number is intrinsic to Keplerian 2D turbulence, and identify it with the 
minimum value for the ratio k^/k y . Note that the latter is qflt max at the time of maximum 
growth, which is limited by viscosity through Eq. (81). The rationale of our argument is that 
if the energy growth for modes with k^/k y = 6 — 11 is significantly hindered by viscosity, 
then vortices of the corresponding scale damp and thus turbulence does not develop at 
smaller scales. Now, defining the effective Reynolds number for a simulation box of size L 
as 1Z = QL 2 /u, and noting that the minimum value of k y available to the modes is 27r/L for 
periodic boundary conditions, Eq. (81) for a 2D Keplerian flow gives 



Based on analogy with the onset of turbulence in subcritical Couette flow, MAN04 find 
a value for 1Z C that is consistent with the high end of the above range. We thus conclude that 
the critical Reynolds number required for sustaining turbulence in 3D numerical simulations 
of a local Keplerian flow is (at least) ~ 10 4 , which is comparable with the effective dynamical 




3 



(99) 
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range in BHSW simulations. Therefore, we expect larger 3D simulations, with appropriately 
chosen initial conditions, and/or long enough run time to start seeing the possible emergence 
of Keplerian turbulence. 

Here, we note that as the turbulent vortices are stretched in the streamwise direction, 
more resolution in the radial (compared to the streamwise) direction is needed to resolve the 
turbulent structure. Therefore, the most efficient way to simulate the turbulent regime is to 
look at a simulation box that is elongated in the streamwise (y) direction by a factor of 6 — 11 
compared to the radial (x) direction, while similar number of (asymmetric) cells in x and 
y directions are used. This will insure that the largest (i.e. least damped) sheared vortices 
can fit in the simulation box, while the vortices are adequately resolved in both directions. 
We do not have a solid answer for the optimum vertical (z) size of the simulation box but 
note that, due to the decay of vorticity in 2D, the vertical structure is necessary to sustain 
the turbulence. 



8. Conclusions 

In this paper, for the first time, we present a three-dimensional analytic study of the 
transient growth of energy for plane wave disturbances of a rotating shear flow, in the limit of 
large speed of sound and small viscosity, and consider its direct application to astrophysical 
accretion disks. We see that, although the growth is fastest for 2D disturbances, a large 
transient energy growth is a generic feature of non-axisymmetric disturbances with relatively 
large radial leading wave numbers (k x ^> k y , k z ; k x k y < 0). The maximum growth is quadratic 
in time for 2D disturbances, but becomes linear in time for relatively large vertical wave 
numbers (k z > k y ), and long times. 

After including the effects of viscosity and compressibility /finite disk thickness, we show 
that the maximum energy growth scales as (1Z/ log72) 2//3 , where, 72. ^> 1, is the Reynolds 
number. Therefore, for neutral astrophysical accretion disks with 1Z ~ 10 10 — 10 14 , assuming 
that there exists a non-linear feedback process to repopulate the growing disturbances, the 
transient growth can act as an alternative to the more conventional MRI instability in ionized 
disks, in starting and sustaining the turbulence necessary to explain the observed accretion 
efficiencies. 

Finally, we present a heuristic argument based on the aspect ratio of sheared turbulent 
vortices to show that the critical Reynolds number for sustaining Keplerian hydrodynamic 
turbulence in a periodic shearing box should be ~ 10 4 . A hydrodynamical simulation needs 
to start with a fine tuned initial condition or last many orbital times to reach the sustained 
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turbulent regime. 

We would like to thank the anonymous referee for helpful comments and suggestions. 
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v v = - qQ x 



Towards Central Object 



Fig. 1. — A plot of the unperturbed flow in the local co moving box studied in this paper. 
The thick arrows represent the velocity field. 
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Fig. 2. — The logarithmic slope of the maximum growth function, a, i.e. G max oc (qflt) a , 
plotted as a function of k z /k y for a Keplerian disk. For larger values of k z /k y , a remains 
constant. 



